// Magneticle, a magnetic field visualization using Biot-Savart law
// Copyright (C) 2005-2008 David J. Beal
// 
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
// 
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
// 
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see <http://www.gnu.org/licenses/>.

#ifdef __USE_ISOC99
#undef __USE_ISOC99
#endif
#define __USE_ISOC99 1
#include <assert.h>
#include <GL/glut.h>
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include <time.h>
#include <sys/types.h>
#include <unistd.h>
#include <vector>
#include <string.h>
using namespace std;

/* might be 4, for example */
#define DIMENSION 3
#define VECFIELD_DIMENSION DIMENSION
#define VFD VECFIELD_DIMENSION

int viewport_width;
int viewport_height;
int draw_stereo;
double stereo_sep;

int animate_curve=1;
int alpha_magnitude;
int alpha_direction=1;
double limd;
int do_perturb;
double perturb;

//double mag_current=0.022;
double mag_current=10.000;
int mag_resolution=16;


typedef struct
{
	//Triplet trip;
	double x[VFD];
	//Triplet x;
	double v[VFD];
	double a[VFD];
	double mass;
	double color[4];
}
point_t;
typedef struct
{
	vector<point_t> points;
}
pointset_t;
pointset_t pointset;


typedef struct
{
	float shininess;
	/* white light component coefficients */
	float specular[4];
	float ambient[4];
	float diffuse[4];
	float emission[4];
	int sides; /* 1 or 2 */
	int orientation; /* +1 or -1 for positive or negative orientation */
	enum
	{
		POINT,
		LINE,
		FILL
	}rasterization_mode;
	enum
	{
		CONSTANT,
		LINEAR,
		QUADRATIC
	}interpolation_model;
}material_t;
material_t materials[1];
typedef struct
{
	float x[4];
	float specular[4];
	float diffuse[4];
}light_t;
light_t lights[1];
int enable_lights;
float ambient_light[4];
typedef struct
{
	/* midpoint rule */
	unsigned int n[VFD];
	double lim[VFD][2];
	void (*f)(double const *,double *,point_t &p);
	double normal_alpha_magnitude;
	double color[4];
	int on;
}vecfield_t;
vecfield_t vecfields[2];


double vector_field_scale;
int draw_vector_field;
int draw_pointset_flag;
int animate_pointset_flag;
int pointset_size=256;
double point_max_radius=128;
//double initial_point_radius=.02;
double initial_point_radius=10.00;
double initial_point_velocity=1;
double const_c;
void draw_point3dv(double const *pvec);
void draw_point3d(double x,double y,double z);
typedef struct
{
	/*
	enum
	{
		Z_FXY=0,
		X_FYZ=1,
		Y_FZX=2,
		PARAM=3,
	}type;
	*/
	/* vector components */
	/* limits */
	double lim[3][2];
	/* number of isolines used to draw surface */
	int n[3];
	/* 
	double (*f)(double const *vec);
	*/
	/* surface function */
	void (*f)(double const *uv,double *r);
	/* unit normal function */
	void (*nf)(double const *uv,double *r);
	/* draw this surface */
	bool draw;
	/* material pointer */
	int material;
}surf_t;
surf_t surfaces[1];

struct Triplet
{
	double x[VFD];
	double &operator[](int offset)
	{
		return this->x[offset];
	}
	operator double const *(void)
	{
		return &this->x[0];
	}
	operator double *(void)
	{
		return &this->x[0];
	}
};

Triplet ppp;
double *thex=&ppp[0];
double *they=ppp;
double const *thez=ppp;

struct Curve;

typedef void (*Curve_function_t)(Curve &c,double,double *);
typedef void (*curve_function_t)(double,double *);

struct Curve
{
	Curve_function_t fn;
	double t0;
	double t1;
	double dt;
	unsigned int nt;
	int enabled;

	double at;
	double adt;

	Curve
	(
		 Curve_function_t _fn
		,double _t0
		,double _t1
		,unsigned int _nt
		,double _at
		,double _adt
	)
		:fn(_fn)
		,t0(_t0)
		,t1(_t1)
		,nt(_nt)
		,enabled(1)
		,at(_at)
		,adt(_adt)
	{
		this->dt=(t1-t0)/nt;
	}

	void biot_savart(double dampen,int n0,int n1,double mt0,double dtm,double const *x,long double *b);

	void draw_curve(void);
};

vector<Curve> curves;

struct MagField
{
	MagField(void)
	{
	}
	vector<Curve> curves;
};

vector<MagField> magfields;

int draw_surface;
double points[1][3]=
{
	{6,4,9}
};
int draw_point;
double plane_scale;
int bool_draw_plane;
int plane_density;
int draw_tan_curve;
int draw_norm_curve;
int draw_param_curve;
int ortho;
double orscale;
double t0;
double t1;
double dt;
double t20;
double t21;
double d2t;
double t30;
double t31;
double d3t;
int curve_resolution;
double tp;
double scale;
double scalex;
double scaley;
double scalez;
double mrotx;
double mroty;
double mrotz;
double mrd;
double eyex;
double eyey;
double eyez;
double cenx;
double ceny;
double cenz;
double upx;
double upy;
double upz;
double fov;
double aspect;
double znear;
double zfar;



/* random long double from 0.0 to 1.0 inclusive */
#define randfloat ((long double)rand()*(1./((long double)RAND_MAX+1.0)))
/* random long double from -1.0 to 1.0 inclusive */
#define randfloatn (randfloat*2.0-1.0)

#define DBUG(VAR) cout<<" ## VAR ## :"<<(VAR)<<endl

#define drawOneLine(x1,y1,x2,y2)  glBegin(GL_LINES);  \
   glVertex2f ((x1),(y1)); glVertex2f ((x2),(y2)); glEnd();

#define ELEMSOF(a) (sizeof(a)/sizeof((a)[0]))
#define CLONE(a,b) typeof a b = a
#define CLONEOF(a,b) CLONE(a,b)

/* C = A x B */
void cross(double const *a,double const *b,double *c)
{
	double x1=a[0];
	double y1=a[1];
	double z1=a[2];
	double x2=b[0];
	double y2=b[1];
	double z2=b[2];
	c[0]=y1*z2-z1*y2;
	c[1]=z1*x2-x1*z2;
	c[2]=x1*y2-y1*x2;
	/*
	c[0]=a[1]*b[2]-a[2]*b[1];
	c[1]=a[2]*b[0]-a[0]*b[2];
	c[2]=a[0]*b[1]-a[1]*b[0];
	*/
}

double dot(double const *a,double const *b)
{
	return a[0]*b[0]+a[1]*b[1]+a[2]*b[2];
}

void normalize(double const *x,double *n)
{
	long double nd=1./sqrt(dot(x,x));
	n[0]=x[0]*nd;
	n[1]=x[1]*nd;
	n[2]=x[2]*nd;
}

void init(void) 
{
	glClearColor(0.0, 0.0, 0.0, 0.0);
}

void reshape(int w,int h)
{
	viewport_width=w;
	viewport_height=h;
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	aspect=(double)viewport_width/viewport_height;
	if(1)
	{
		fprintf(stderr,"reshape viewport to width %d height %d aspect %f\n",w,h,aspect);
	}
}

/* depends on curve_resolution */
void calc_t()
{
	t0=0;
	t1=2*M_PI;
	dt=(t1-t0)/curve_resolution;
	tp=0;

	t20=0;
	t21=2*M_PI;
	d2t=(t21-t20)/curve_resolution;

	t30=0;
	t31=2*M_PI;
	d3t=(t31-t30)/curve_resolution;
}

#define SETV(V,X,Y,Z) (V)[0]=(X);(V)[1]=(Y);(V)[2]=(Z)
#define SET4(V,X,Y,Z,W) (V)[0]=(X);(V)[1]=(Y);(V)[2]=(Z);(V)[3]=(W)
double G=6.673e-11;

void nullfield(double const *x,double *v)
{
	SETV(v,0,0,0);
}

//double grav_threshold=1.125;
double grav_threshold=0.200;
double universe_dampen=0.990;
//double universe_dampen=1.000;
double grav_center[3]={0,0,0};
double grav_mass=100000000000.0;

void remove_point(pointset_t &pointset,unsigned int index);
void initialize_point(point_t &point);
void gravfield(double const *ix,double *v,point_t &point)
{
	/* domain set of points */
	double x[3];
	int n;
	n=VFD;
	double d2=0;
	while(n--)
	{
		x[n]=grav_center[n]-ix[n];
		d2+=x[n]*x[n];
	}
	double d2_32=pow(d2,3./2);
	double factor=grav_mass*G/d2_32;
	if(1)
	{
		if(d2<grav_threshold)
		{
			#if 0
			if(0)
			{
				if(0)
				{
				}
				else
				{
					if(0)
					{
						point.v[0]*=-1;
						point.v[1]*=-1;
						point.v[2]*=-1;
					}
				}
				initialize_point(point);
			}
			remove_point(point);
			#endif
		}
		else
		{
			n=VFD;
			double b[3];
			while(n--)
			{
				b[n]=x[n]*factor;
				v[n]=b[n];
			}
		}
	}
}

void draw_vecfield(vecfield_t &vf)
{
	point_t point;
	#if 0
	double x[VFD];
	#endif
	double *x=point.x;
	double v[VFD];
	double r[VFD];
	double y[VFD];
	double b[VFD];
	double d[VFD];
	int m[VFD];
	int l;

	if(perturb)
	{
		srand(getpid());
	}

	l=VFD;
	while(l--)
	{
		d[l]=(vf.lim[l][1]-vf.lim[l][0])/vf.n[l];
		b[l]=vf.lim[l][0]+d[l]*.5;
		m[l]=vf.n[l];
	}

	glEnable(GL_LINE_SMOOTH);
	glEnable(GL_BLEND);
	glBlendFunc(GL_SRC_ALPHA,GL_ONE);

	glBegin(GL_LINES);

	double vfc[4]={vf.color[0],vf.color[1],vf.color[2],vf.color[3]};

	m[2]=vf.n[2];
	while(m[2]--)
	{
		m[1]=vf.n[1];
		while(m[1]--)
		{
			m[0]=vf.n[0];
			while(m[0]--)
			{
				l=VFD;
				while(l--)
				{
					x[l]=b[l]+m[l]*d[l];
					if(perturb)
					{
						x[l]+=d[l]*perturb*randfloatn*randfloatn;
					}
				}

				point_t point;
				vf.f(x,v,point);

				if(alpha_magnitude)
				{
					double m=(dot(v,v));
					glColor4d(vf.color[0],vf.color[1],vf.color[2],fmin(m/vf.normal_alpha_magnitude,1));
					m=1./sqrt(m);
					v[0]*=m;
					v[1]*=m;
					v[2]*=m;
				}
				else
				{
					normalize(v,v);
				}

				v[0]*=vector_field_scale;
				v[1]*=vector_field_scale;
				v[2]*=vector_field_scale;

				y[0]=x[0]-v[0]*.5;
				r[0]=x[0]+v[0]*.5;
				y[1]=x[1]-v[1]*.5;
				r[1]=x[1]+v[1]*.5;
				y[2]=x[2]-v[2]*.5;
				r[2]=x[2]+v[2]*.5;

				if(1)
				{
					glColor4dv(vfc);
				}
				if(alpha_direction)
				{
				}
				glVertex3dv(y);
				if(alpha_direction)
				{
					glColor4d(0,0,0,0);
				}
				glVertex3dv(r);
			}
		}
	}
	glEnd();

	if(perturb)
	{
		static int x=0;
		srand(getpid()+x++);
	}
}

void draw_pointset(pointset_t const &pointset)
{
	int num_points=pointset.points.size();
	int n=num_points;
	CLONEOF(&pointset.points[0],pvec);

	glEnable(GL_LINE_SMOOTH);
	glEnable(GL_BLEND);
	glBlendFunc(GL_SRC_ALPHA,GL_ONE);

	glBegin(GL_LINES);
	while(n--)
	{
		CLONEOF(pvec[n],&pv);
		double const *x=pv.x;
		double const *v=pv.v;
		//double const *a=pv.a;

		glColor4dv(pv.color);
		glVertex3dv(x);
		glVertex3d(x[0]+v[0],x[1]+v[1],x[2]+v[2]);
	}
	glEnd();

	if(1)
	{
		glEnable(GL_POINT_SMOOTH);
		glBegin(GL_POINTS);
		n=num_points;
		while(n--)
		{
			CLONEOF(pvec[n],&pv);
			double const *x=pv.x;
			double const *v=pv.v;

			double r2=0;
			r2+=v[0]*v[0];
			r2+=v[1]*v[1];
			r2+=v[2]*v[2];

			if(r2<1)
			{
				glColor4dv(pv.color);
				glVertex3dv(x);
			}
		}
		glEnd();
	}
}
void animate_pointset_through_vecfield(vecfield_t &vf,pointset_t &pointset)
{
	int num_points=pointset.points.size();
	int n=num_points;
	CLONEOF(&pointset.points[0],pvec);

	while(n--)
	{
		CLONEOF(pvec[n],&pv);
		//double *pv=pvec[n];
		double *x=pv.x;
		double *v=pv.v;
		double *a=pv.a;
		double field_accel_out[3];

		vf.f(x,&field_accel_out[0],pv);

		if(1)
		{
			int n=VFD;
			while(n--)
			{
				field_accel_out[n]*=vector_field_scale;
				a[n]+=field_accel_out[n];
				//v[n]+=a[n];
				//v[n]*=universe_dampen;
				//x[n]+=v[n];
			}
		}
	}
}
void apply_acceleration(pointset_t &pointset)
{
	int num_points=pointset.points.size();
	int n=num_points;
	CLONEOF(&pointset.points[0],pvec);

	while(n--)
	{
		CLONEOF(pvec[n],&pv);
		//double *pv=pvec[n];
		double *x=pv.x;
		double *v=pv.v;
		double *a=pv.a;

		if(1)
		{
			int n=VFD;
			while(n--)
			{
				//a[n]*=universe_dampen;
				v[n]+=a[n];
				a[n]=0;
				//v[n]*=universe_dampen;
				x[n]+=v[n];
				v[n]=0;
			}
		}
	}
}
void stop_pointset(pointset_t &pointset)
{
	int num_points=pointset.points.size();
	while(num_points--)
	{
		int i=3;
		while(i--)
		{
			pointset.points[num_points].v[i]=0;
		}
	}
}
void clear_pointset_acceleration(pointset_t &pointset)
{
	int num_points=pointset.points.size();
	while(num_points--)
	{
		int n=VFD;
		while(n--)
		{
			pointset.points[num_points].a[n]=0;
		}
	}
}


void paramfunc1(double t,double *pvec)
{
	// SETV(pvec,0.8*cos(t),0.8*sin(t),0.8*t/t0);
	// SETV(pvec,0.8*cos(t),0.8*sin(t),0.8*cos(t*3));
	// SETV(pvec,0,0,t);
	// SETV(pvec,0.8*cos(t),0.8*sin(t),t/t0);
	//
	SETV(pvec,0.8*cos(t),0.8*sin(t),0);
}
void paramfunc2(double t,double *pvec)
{
	SETV(pvec,0.8*cos(t),0,0.8*sin(t));
}
void trefoil_knot(double t,double *pvec)
{
	SETV(pvec,(2+cos(1.5*t))*cos(t),(2+cos(1.5*t))*sin(t),sin(1.5*t));
}
void anti_circle(Curve &c,double t,double *pvec)
{
	SETV(pvec,2*cos(t),-4+2*sin(c.at),2*sin(t));
}
void circle(Curve &c,double t,double *pvec)
{
	SETV(pvec,cos(-t),-1,sin(-t));
}
void circle0(Curve &c,double t,double *pvec)
{
	SETV(pvec,cos(t),1,sin(t));
}
void circle0_anim(Curve &c,double t,double *pvec)
{
	double at=c.at;
	SETV(pvec,sin(at)*cos(t),1*cos(at),sin(t));
}
void circle1(Curve &c,double t,double *pvec)
{
	SETV(pvec,cos(c.at)*cos(t),sin(c.at)*sin(t),1);
}
void circle2(Curve &c,double t,double *pvec)
{
	SETV(pvec,cos(c.at)*cos(-t),sin(c.at)*sin(-t),-1);
}
void unittan1(double t,double *pvec)
{
	SETV(pvec,1,1/t,8*t/pow(pow(t,2)+1,2));
}
void unitnorm1(double t,double *pvec)
{
	SETV(pvec,0,-1/pow(t,2),8*(1-3*pow(t,2))/pow(pow(t,2)+1,3));
}

void spiral_0(Curve &c,double t,double *pvec)
{
	double rotations=4;
	double const radius=2*t;
	double theta=-t*2*M_PI*rotations;
	double x=radius*cos(theta);
	double y=4-t;
	double z=radius*sin(theta);
	double ox=x;
	double oy=y;
	if(animate_curve)
	{
		x=ox*cos(c.at)+oy*sin(c.at);
		y=ox*sin(c.at)-oy*cos(c.at);
	}
	SETV(pvec,x,y,z);
}

void spiral_1(Curve &c,double t,double *pvec)
{
	double rotations=4;
	double const radius=2*t;
	double theta=-t*2*M_PI*rotations;
	double x=radius*cos(-theta);
	double y=-4+t;
	double z=radius*sin(-theta);
	double ox=x;
	double oy=y;
	if(animate_curve)
	{
		x=ox*cos(c.at)+oy*sin(c.at);
		y=ox*sin(c.at)-oy*cos(c.at);
	}
	SETV(pvec,x,y,z);
}

void toroidal_spiral(Curve &c,double t,double *pvec)
{
	#if 0
	double rotations=4;
	double const radius=2*t;
	double theta=-t*2*M_PI*rotations;
	double x=radius*cos(-theta);
	double y=-4+t;
	double z=radius*sin(-theta);
	double ox=x;
	double oy=y;
	if(animate_curve)
	{
		x=ox*cos(c.at)+oy*sin(c.at);
		y=ox*sin(c.at)-oy*cos(c.at);
	}
	#endif
	double x=0;
	double y=0;
	double z=0;
	double const ts=20;
	double const s=8;
	x=(4+sin(ts*t))*cos(t);
	y=(4+sin(ts*t))*sin(t);
	z=cos(ts*t);
	x*=s;
	y*=s;
	z*=s;
	x+=s*4;
	//y+=s*2;
	//z+=s*2;
	// trying to place the toroidal spiral so that the particles are inside the tube and begin cycling
	SETV(pvec,x,y,z);
}

void cardioid_0(Curve &c,double t,double *pvec)
{
	double const radius=2*sin(t*M_PI);
	double theta=t*2*M_PI;
	double x=radius*cos(theta);
	double y=0;
	double z=radius*sin(theta);
	SETV(pvec,x,y,z);
}

void initialize_curves()
{
	curves.clear();
	curve_resolution=mag_resolution;
	double cadt=0.08;
	if(0)
	{
		curves.push_back(Curve(circle,0,2*M_PI,curve_resolution,0,cadt));
	}
	if(0)
	{
		curves.push_back(Curve(anti_circle,0,2*M_PI,curve_resolution,0,cadt));
	}
	if(0)
	{
		curves.push_back(Curve(circle0_anim,0,2*M_PI,curve_resolution,0,cadt));
	}
	if(0)
	{
		curves.push_back(Curve(circle1,0,2*M_PI,curve_resolution,0,cadt));
	}
	if(0)
	{
		curves.push_back(Curve(circle2,0,2*M_PI,curve_resolution,0,cadt));
	}
	if(0)
	{
		curves.push_back(Curve(spiral_0,0,1,curve_resolution,0,cadt));
	}
	if(0)
	{
		curves.push_back(Curve(spiral_1,0,1,curve_resolution,0,cadt));
	}
	if(1)
	{
		curves.push_back(Curve(toroidal_spiral,0,2*M_PI,curve_resolution,0,cadt));
	}
	//curves.push_back(Curve(trefoil_knot,0,4*M_PI,curve_resolution,0,cadt));
}


void Curve::biot_savart(double dampen,int n0,int n1,double mt0,double dtm,double const *x,long double *b)
{
	double r[3];
	double p[3];
	double dr[3];
	double pn[3];
	double dc[3];
	double d;

	this->fn(*this,mt0+(n0)*dtm,r);
	this->fn(*this,mt0+(n1)*dtm,dr);

	dr[0]-=r[0];
	dr[1]-=r[1];
	dr[2]-=r[2];

	p[0]=x[0]-r[0];
	p[1]=x[1]-r[1];
	p[2]=x[2]-r[2];

	normalize(p,pn);
	cross(dr,pn,dc);

	d=dtm/dot(p,p);

	d*=mag_current*dampen;

	b[0]+=dc[0]*d;
	b[1]+=dc[1]*d;
	b[2]+=dc[2]*d;
}

#if 0
void biot_savart(void (*pf)(double,double *),double dampen,int n0,int n1,double mt0,double dtm,double const *x,long double *b)
{
	double r[3];
	double p[3];
	double dr[3];
	double pn[3];
	double dc[3];
	double d;

	pf(mt0+(n0)*dtm,r);
	pf(mt0+(n1)*dtm,dr);

	dr[0]-=r[0];
	dr[1]-=r[1];
	dr[2]-=r[2];

	p[0]=x[0]-r[0];
	p[1]=x[1]-r[1];
	p[2]=x[2]-r[2];

	normalize(p,pn);
	cross(dr,pn,dc);

	d=dtm/dot(p,p);

	d*=mag_current*dampen;

	b[0]+=dc[0]*d;
	b[1]+=dc[1]*d;
	b[2]+=dc[2]*d;
}
#endif

#if 0
/* assumes paramfunc1 */
void magfield1(double const *x,double *v)
{
	long double b[3];
	/*
	*/

	/* uses Biot-Savart Law */
	/* obtain line integral for magnetic field at point given by argument vector x */

	int n;
	double dtm;
	double d3tm;
	double idt;

	b[0]=0;
	b[1]=0;
	b[2]=0;

	if(0)
	{
		n=abs((int)(const_c));
		//n=1;
	}
	else
	{
		n=mag_resolution;
	}

	idt=1./n;

	dtm=(t1-t0)/n;
	d3tm=(t31-t30)/n;
	while(n--)
	{
		//biot_savart(paramfunc1,idt,n+1,n,t0,dtm,x,b);
		//biot_savart(paramfunc2,idt,n+1,n,t0,dtm,x,b);
		//biot_savart(trefoil_knot,idt,n+1,n,t30,d3t,x,b);
		biot_savart(circle,idt,n+1,n,t30,d3tm,x,b);
	}

	SETV(v,b[0],b[1],b[2]);
}
#endif

#if 0
/* assumes paramfunc1 */
void magfield2(double const *x,double *v)
{
	long double b[3];

	int n;
	double dtm;
	double d3tm;
	double idt;

	b[0]=0;
	b[1]=0;
	b[2]=0;

	if(0)
	{
		n=abs((int)(const_c));
		//n=1;
	}
	else
	{
		n=mag_resolution;
	}

	idt=1./n;

	dtm=(t1-t0)/n;
	d3tm=(t31-t30)/n;
	while(n--)
	{
		//biot_savart(paramfunc1,idt,n+1,n,t0,dtm,x,b);
		//biot_savart(paramfunc2,idt,n+1,n,t0,dtm,x,b);
		//biot_savart(trefoil_knot,idt,n+1,n,t30,d3t,x,b);
		biot_savart(circle0,idt,n+1,n,0,2*M_PI*idt,x,b);
	}

	SETV(v,b[0],b[1],b[2]);
}
#endif

void magfield3(double const *x,double *v,point_t &point)
{
	long double b[]={0,0,0};
	unsigned int n=mag_resolution;
	double idt=1./n;

	while(n--)
	{
		//biot_savart(circle0,idt,n+1,n,0,2*M_PI*idt,x,b);
		int c=curves.size();
		while(c--)
		{
			Curve &curve=curves[c];
			if(curve.enabled)
			{
				double dtm=(curve.t1-curve.t0)*idt;
				curve.biot_savart(idt,n+1,n,curve.t0,curve.t1*idt,x,b);
			}
		}
	}

	SETV(v,b[0],b[1],b[2]);
}

double planefunc1xy(double const *x)
{
	/* cx + y + z = 1 */
	/* z = 1 - cx - y */
	return 1-const_c*x[0]-x[1];
}

double planefunc2xy(double const *x)
{
	/* x - cy + cz = -1 */
	/* z = (cy - x - 1)/c */
	return (const_c*x[1] - x[0] - 1)/const_c;
}

double xyplane(double const *x)
{
	return 0;
}

double surffunc1xfyz(double const *x)
{
	double const sdevx=0.15;
	double const sdevy=0.15;
	double const f=1./(2*M_PI*sdevx*sdevy);
	double const meanx=1*const_c;
	double const meany=1;
	return f*exp(-pow(x[0]-meanx,2)-pow(x[1]-meany,2));
}

double surffunc2zfxy(double const *x)
{
	return 2-x[1];
}

double surf4(double const *x)
{
	return 2-pow(x[0],4)-pow(x[1],4);
}

/* u[2] */
/* r=f(u0,u1) */
void paramsurf(double const *u,double *r)
{
	SETV(r,u[0],u[1],2-pow(u[0],4)-pow(u[1],4));
}
void paramnorm(double const *u,double *n)
{
	SETV(n,4*pow(u[0],3),4*pow(u[1],3),1);
	double nf=1./sqrt(n[0]*n[0]+n[1]*n[1]+n[2]*n[2]);
	n[0]*=nf;
	n[1]*=nf;
	n[2]*=nf;
}

void paramsphere(double const *u,double *x)
{
	SETV(x,const_c*sin(u[0])*cos(u[1]),const_c*sin(u[0])*sin(u[1]),const_c*cos(u[0]));
	SETV(x,sin(u[0])*cos(u[1]),cos(2*u[1])*sin(3*u[0]),cos(u[1]));
}

void spherenorm(double const *u,double *x)
{
	SETV(x,sin(u[0])*cos(u[1]),sin(u[0])*sin(u[1]),cos(u[0]));
}

void prepare_surface(surf_t &s)
{
#if 1
	CLONE(materials[s.material],&mat);
	assert(mat.sides==1 || mat.sides==2);
	if(enable_lights)
	{
		if(mat.sides==1)
		{
			glLightModelf(GL_LIGHT_MODEL_TWO_SIDE,0);
		}
		else
		{
			glLightModelf(GL_LIGHT_MODEL_TWO_SIDE,1);
		}

		if(1)
		{
			glEnable(GL_LIGHTING);
			glEnable(GL_LIGHT0);
			glLightModelfv(GL_LIGHT_MODEL_AMBIENT,ambient_light);
			glLightModelf(GL_LIGHT_MODEL_LOCAL_VIEWER,0);

			glLightfv(GL_LIGHT0,GL_DIFFUSE,lights[0].diffuse);
			glLightfv(GL_LIGHT0,GL_SPECULAR,lights[0].specular);

			glEnable(GL_AUTO_NORMAL);
			glEnable(GL_NORMALIZE);
		}



		if(1)
		{
			GLenum e=GL_FRONT_AND_BACK;
			glMaterialfv(e,GL_AMBIENT,mat.ambient);
			glMaterialfv(e,GL_DIFFUSE,mat.diffuse);
			glMaterialfv(e,GL_SPECULAR,mat.specular);
			glMaterialfv(e,GL_EMISSION,mat.emission);
			glMaterialf(e,GL_SHININESS,mat.shininess);
			/*
			glColorMaterial(e,GL_DIFFUSE);
			glEnable(GL_COLOR_MATERIAL);
			*/
		}
	}
	switch(mat.interpolation_model)
	{
	    case material_t::CONSTANT:
		glShadeModel(GL_FLAT);
		break;
	    case material_t::LINEAR:
		glShadeModel(GL_SMOOTH);
		break;
	    case material_t::QUADRATIC:
		/* ?? */
		break;
	}
	switch(mat.rasterization_mode)
	{
	    case material_t::POINT:
		glPolygonMode(GL_FRONT_AND_BACK,GL_POINT);
		break;
	    case material_t::LINE:
		/* not doing polygons */
		glPolygonMode(GL_FRONT_AND_BACK,GL_LINE);
		glDisable(GL_DEPTH_TEST);
		glEnable(GL_BLEND);
		glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA);
		glEnable(GL_LINE_SMOOTH);
		glLineWidth(1.0);
		break;
	    case material_t::FILL:
		glPolygonMode(GL_FRONT_AND_BACK,GL_FILL);
		glEnable(GL_DEPTH_TEST);
		glDepthFunc(GL_LESS);
		glDisable(GL_BLEND);
		break;
	}
#endif
}

void deprepare_surface(surf_t &s)
{
	if(enable_lights)
	{
		glDisable(GL_LIGHTING);
		glDisable(GL_AUTO_NORMAL);
		glDisable(GL_NORMALIZE);
	}
	glPolygonMode(GL_FRONT_AND_BACK,GL_LINE);
	glDisable(GL_DEPTH_TEST);
	glEnable(GL_BLEND);
	glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA);
	glEnable(GL_LINE_SMOOTH);
	glLineWidth(1.0);
}

void rasterize_surface_lines(surf_t &s)
{
	int i,j,k;
	double b[2];
	double d[2];
	double u[2];
	double x[3];
	int const lo=1; /* line offset */

	k=-1;
	while(++k<2)
	{
		d[k]=(s.lim[k][1]-s.lim[k][0])/s.n[k];
		b[k]=s.lim[k][0];
	}

	k=-1;
	while(++k<2)
	{
		/* number of v */
		j=-1;
		while(++j<s.n[k^1]+lo)
		{
			u[k^1]=b[k^1]+j*d[k^1];
			glBegin(GL_LINE_STRIP);
			/* number of u */
			i=-1;
			while(++i<s.n[k^0]+lo)
			{
				u[k^0]=b[k^0]+i*d[k^0];
				s.f(u,x);
				glVertex3dv(x);
			}
			glEnd();
		}
	}
}

void rasterize_surface_quads(surf_t &s)
{
	int i,j,k;
	double b[2];
	double d[2];
	double u[2];
	double x[3];
	int ND=1;

	k=-1;
	while(++k<2)
	{
		d[k]=(s.lim[k][1]-s.lim[k][0])/s.n[k];
		b[k]=s.lim[k][0];
	}

	k=-1;
	while(++k<ND)
	{
		/* number of v */
		j=-1;
		while(++j<s.n[k^1])
		{
			glBegin(GL_QUAD_STRIP);
			/* number of u */
			i=-1;
			while(++i<s.n[k^0]+1)
			{
				u[k^0]=b[k^0]+i*d[k^0];
				u[k^1]=b[k^1]+j*d[k^1];
				s.f(u,x);
				if(1)
				{
					//s.nf(u,n);
					//glNormal3dv(n);
				}
				glVertex3dv(x);
				if(1)
				{
					u[k^1]=b[k^1]+(j+1)*d[k^1];
					//s.nf(u,n);
					//glNormal3dv(n);
					s.f(u,x);
					glVertex3dv(x);
				}
			}
			glEnd();
		}
	}

#if 0
	double v[3];
	int i,j,k;
	double d[3];

	k=3;
	while(k--)
	{
		d[k]=(s.x[k][1]-s.x[k][0])/s.n[k];
	}

	/* u and v */
	k=2;
	while(k--)
	{
		/* c=f(a,b) */
		int q=s.type;
		int a=(q+(k^1))%3;
		int b=(q+(k^0))%3;
		int c=(q+2)%3;
		j=s.n[a];
		while(j--)
		{
			/* draw y=k traces */
			glBegin(GL_QUADS);
			v[a]=s.x[a][0]+j*d[a];
			i=s.n[b];
			while(i--)
			{
				v[b]=s.x[b][0]+i*d[b];
				v[c]=s.f(v);
				if(isnormal(v[c]) || v[c]==0.0)
				{
					glVertex3dv(v);
				}
				else
				{
					glEnd();
					glBegin(GL_LINE_STRIP);
				}
			}

			glEnd();
		}
	}
#endif
}

/* draw only parametric surfaces */
void draw_surf(surf_t &s)
{
	/* ends up computing the function 10 billion times more than necessary
	 * why doesn't GL have a surface ability? */

	prepare_surface(s);

	switch(materials[s.material].rasterization_mode)
	{
	    case material_t::FILL:
		rasterize_surface_quads(s);
		break;
	    case material_t::POINT:
	    case material_t::LINE:
		rasterize_surface_lines(s);
		break;
	}

	deprepare_surface(s);
}

double arclength(void (*pf)(double,double *),double t0,double t1,double dt)
{
	double ta;
	double tb;
	long double L=0;
	double pvec0[3];
	double pvec1[3];
	int i;
	int n;

	i=0;
	n=curve_resolution;
	while(n--)
	{
		ta=t0+(i+0)*dt;
		tb=t0+(i+1)*dt;
		pf(ta,&pvec0[0]);
		pf(tb,&pvec1[0]);
		L+=sqrt(pow(pvec1[0]-pvec0[0],2)+
			pow(pvec1[1]-pvec0[1],2)+
			pow(pvec1[2]-pvec0[2],2));
		++i;
	}
	return L;
}

void show_arclength(void (*pf)(double,double *),double t0,double t1,double dt)
{
	double L=arclength(pf,t0,t1,dt);
	printf("arclength[n=%d]=%.15f\n",curve_resolution,L);
}

void Curve::draw_curve(void)
{
		//draw_curve(curve.fn,curve.t0,curve.t1,curve.dt);
		//
	double t0=this->t0;
	double t1=this->t1;
	double dt=this->dt;
	double t=t0;
	double pvec[3];
	int i;
	int n;

	glColor3f(1,1,1);
	glBegin(GL_LINE_STRIP);

	i=0;
	n=curve_resolution+1;
	while(n--)
	{
		t=t0+i*dt;
		this->fn(*this,t,&pvec[0]);
		glVertex3dv(pvec);
		++i;
	}

	glEnd();

	if(animate_pointset_flag)
	{
		this->at+=this->adt;
	}
}
void draw_curve(curve_function_t pf,double t0,double t1,double dt)
{
	double t=t0;
	double pvec[3];
	int i;
	int n;

	glColor3f(1,1,1);
	glBegin(GL_LINE_STRIP);

	i=0;
	n=curve_resolution+1;
	while(n--)
	{
		t=t0+i*dt;
		pf(t,&pvec[0]);
		glVertex3dv(pvec);
		++i;
	}

	glEnd();
}
#if 0
#endif

void draw_curve_object(Curve &curve)
{
	if(curve.enabled)
	{
		curve.draw_curve();
		//draw_curve(curve.fn,curve.t0,curve.t1,curve.dt);
	}
}


void draw_unit_vector(void (*pf)(double,double *),
		      void (*uf)(double,double *),double t)
{
	double pvec[3];
	pf(t,pvec);
	glPushMatrix();
	glTranslated(pvec[0],pvec[1],pvec[2]);

	glBegin(GL_LINES);
	uf(t,pvec);
	glVertex3f(0,0,0);
	glVertex3dv(&pvec[0]);
	glEnd();

	glPopMatrix();
}

void draw_point3dv(double const *pvec)
{
	glPushMatrix();
	glTranslated(pvec[0],pvec[1],pvec[2]);
	glScalef(0.10,0.10,0.10);
	glutWireIcosahedron();
	glPopMatrix();
}

void draw_point3d(double x,double y,double z)
{
	double pvec[3]={x,y,z};
	draw_point3dv(pvec);
}

void draw_the_point(void (*pf)(double,double *),double t)
{
	double pvec[3];
	pf(t,pvec);
	draw_point3dv(pvec);
}

/* draws a random dot that lies on the plane */
void planefunc1(double *pvec)
{
	double z=randfloatn*plane_scale;
	double y=randfloatn*plane_scale;
	double x=6*M_PI-6*y;
	SETV(pvec,x,y,z);
}
void planefunc2(double *pvec)
{
	double z=randfloatn*plane_scale;
	double x=randfloatn*plane_scale;
	double y=M_PI+6*x;
	SETV(pvec,x,y,z);
}
void draw_plane(void (*pf)(double *))
{
	int n=plane_density;
	double pvec[3];
	glBegin(GL_POINTS);
	while(n--)
	{
		pf(&pvec[0]);
		glVertex3dv(pvec);
	}
	glEnd();
}

void draw_axes()
{
	/* draw axes */
	glBegin(GL_LINES);
	glColor3f(1,0,0);
	glVertex3f(0,0,0);
	glVertex3f(1,0,0);
	glColor3f(0,1,0);
	glVertex3f(0,0,0);
	glVertex3f(0,1,0);
	glColor3f(0,0,1);
	glVertex3f(0,0,0);
	glVertex3f(0,0,1);
	glEnd();
}

void draw_unit_circle()
{
	int ns=12;
	int n=ns;
	glBegin(GL_LINE_STRIP);
	while(n--)
	{
		double t=n*(1./ns);
		glVertex3d(cos(t*M_PI*2),sin(t*M_PI*2),0);
	}
	glEnd();
}

void seed_random()
{
	srand(getpid()+time(0));
}

void display(void)
{
	int eye_index=1;

	glClear(GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT);

	if(draw_stereo)
	{
		eye_index=1;
	}

stereo_begin:

	glMatrixMode(GL_PROJECTION);
	glPushMatrix();


	if(ortho)
	{
		glOrtho(-orscale,orscale,-orscale,orscale,-zfar,zfar);
	}
	else
	{
		gluPerspective(fov,aspect*(1./(draw_stereo+1)),znear,zfar);
		//glTranslatef(-eyex,-eyey,-eyez);
		if(draw_stereo)
		{
			double s=-stereo_sep+2*stereo_sep*eye_index;
			glRotatef(s,0,1,0);
		}
		gluLookAt(eyex,eyey,eyez,cenx,ceny,cenz,upx,upy,upz);
	}

	if(draw_stereo)
	{
		glViewport(eye_index*viewport_width/2, 0,viewport_width/2,viewport_height);
	}
	else
	{
		glViewport(0,0,viewport_width,viewport_height);
	}

	glMatrixMode(GL_MODELVIEW);
	glPushMatrix();

	if(enable_lights)
	{
		glLightfv(GL_LIGHT0,GL_POSITION,lights[0].x);
	}

	glRotatef(mrotx,1,0,0);
	glRotatef(mroty,0,1,0);
	glRotatef(mrotz,0,0,1);
	glScalef(scale*scalex,scale*scaley,scale*scalez);


	/* draw axes */
	draw_axes();


	if(draw_param_curve)
	{
		glColor3f(0.5,0.5,0.5);
		//draw_curve(paramfunc1,t0,t1,dt);
		//draw_curve(paramfunc2,t20,t21,d2t);
		//draw_curve(trefoil_knot,t30,t31,d3t);
		//draw_curve(circle,t30,t31,d3t);
		//draw_curve(circle0,t30,t31,d3t);
		int num_curves=curves.size();
		while(num_curves--)
		{
			Curve &curve=curves[num_curves];
			if(curve.enabled)
			{
				curve.draw_curve();
			}
			//draw_curve_object(curves[num_curves]);
		}
	}
	/* draw some given points */
	if(draw_point)
	{
		int n=ELEMSOF(points);
		while(n--)
		{
			draw_point3dv(points[n]);
		}
	}

	if(draw_surface)
	{
		int n=ELEMSOF(surfaces);
		glColor3f(1,1,1);
		while(n--)
		{
			draw_surf(surfaces[n]);
		}
	}

	#if 0
	#endif
	if(draw_tan_curve)
	{
		glColor3f(0,1,1);
		draw_curve(unittan1,t0,t1,dt);
		draw_the_point(unittan1,tp);
	}
	if(draw_norm_curve)
	{
		glColor3f(1,1,0);
		draw_curve(unitnorm1,t0,t1,dt);
		draw_the_point(unitnorm1,tp);
	}

	if(bool_draw_plane)
	{
		glColor3f(1,.5,1);
		draw_plane(planefunc1);
		draw_plane(planefunc2);
	}

	if(draw_vector_field)
	{
		int l=ELEMSOF(vecfields);
		while(l--)
		{
			if(vecfields[l].on)
			{
				draw_vecfield(vecfields[l]);
			}
		}
	}
	if(draw_pointset_flag)
	{
		if(animate_pointset_flag)
		{
			if(1)
			{
				clear_pointset_acceleration(pointset);
			}
			int l=ELEMSOF(vecfields);
			while(l--)
			{
				if(vecfields[l].on)
				{
					animate_pointset_through_vecfield(vecfields[l],pointset);
				}
			}
			apply_acceleration(pointset);
			void reset_out_of_bounds(pointset_t &pointset);
			reset_out_of_bounds(pointset);
			animate_pointset_flag=0;
		}
		if(draw_pointset_flag)
		{
			int l=ELEMSOF(vecfields);
			while(l--)
			{
				draw_pointset(pointset);
			}
		}
	}


	glPopMatrix();
	glMatrixMode(GL_PROJECTION);
	glPopMatrix();
	glMatrixMode(GL_MODELVIEW);

	if(draw_stereo)
	{
		if(eye_index--)
		{
			goto stereo_begin;
		}
	}

	glFlush();
	glutSwapBuffers();
}

void define_surface()
{
	CLONE(surfaces[0],&s0);
	s0.f=paramsphere;
	s0.nf=spherenorm;
	s0.lim[0][0]=0;
	s0.lim[0][1]=M_PI;
	s0.lim[1][0]=0;
	s0.lim[1][1]=2*M_PI;
	s0.n[0]=16;
	s0.n[1]=16;
	s0.n[2]=16;
	s0.material=0;
}

void define_vector_field()
{
	int n=ELEMSOF(vecfields);

	while(n--)
	{
		vecfields[n].on=1;
		vecfields[n].lim[0][0]=-8;
		vecfields[n].lim[0][1]=8;
		vecfields[n].lim[1][0]=-8;
		vecfields[n].lim[1][1]=8;
		vecfields[n].lim[2][0]=-8;
		vecfields[n].lim[2][1]=8;
		vecfields[n].n[0]=8;
		vecfields[n].n[1]=8;
		vecfields[n].n[2]=8;
		vecfields[n].normal_alpha_magnitude=1;
		if(0)
		{
			vecfields[n].color[0]=0.8;
			vecfields[n].color[1]=.243;
			vecfields[n].color[2]=.421;
			vecfields[n].color[3]=.625;
		}
		else
		{
			vecfields[n].color[0]=.101;
			vecfields[n].color[1]=.826;
			vecfields[n].color[2]=.214;
			vecfields[n].color[3]=.625;
		}
		normalize(vecfields[n].color,vecfields[n].color);
	}

	vecfields[0].f=magfield3;
	vecfields[1].f=gravfield;

	if(1)
	{
		vecfields[1].on=0;
	}
}

void initialize_point(point_t &point)
{
	#if 0
	CLONEOF(pointset.points[num_points],&point);
	#endif
	int n=VFD;
	double r=initial_point_radius*randfloat;
	double t=randfloat*M_PI*2;
	double p=randfloat*M_PI*2;
	point.x[0]=r*sin(p)*cos(t);
	point.x[1]=r*sin(p)*sin(t);
	point.x[2]=r*cos(p);
	while(n--)
	{
		//point.x[n]=0;
		//point.v[n]=0;
		point.a[n]=0;
		//point.v[n]=0;
		point.v[n]=randfloatn*randfloat*randfloat*initial_point_velocity;
		//point.v[n]=cos(randfloat*M_PI*2)*randfloat*.01;
		//point.a[n]=randfloatn*randfloatn*randfloatn*.01;
	}
	n=4;
	while(n--)
	{
		point.color[n]=fmin(0.75+randfloat*.25,1);
		//point.color[n]=1.0;
	}
	point.color[3]=0.6;
}
void add_points(pointset_t &pointset,int count)
{
	while(count--)
	{
		point_t point;
		pointset.points.push_back(point);
		initialize_point(pointset.points.back());
	}
}
void remove_point(pointset_t &pointset,unsigned int index)
{
	if(index<pointset.points.size())
	{
		vector<point_t>::iterator it=pointset.points.begin()+index;
		pointset.points.erase(it);
	}
}
void remove_points(pointset_t &pointset,int count)
{
	while(count--)
	{
		if(pointset.points.size())
		{
			remove_point(pointset,rand()%pointset.points.size());
		}
		else
		{
			count=0;
		}
	}
	#if 0
	vector<point_t>::iterator it=pointset.points.begin();
	int l=pointset.points.size();
	if(l)
	{
		while(count)
		{
			if(rand()&1)
			{
				vector<point_t>::iterator nit=it;
				pointset.points.erase(it);
				it=nit++;
				count-=1;
			}
			else
			{
				it++;
			}
			if(it==pointset.points.end())
			{
				it=pointset.points.begin();
			}
		}
	}
	#endif
	#if 0
	while(count--)
	{
		if(pointset.points.size())
		{
			#if 0
			if(rand()&1)
			{
			}
			int rp=rand()%pointset.points.size();
			pointset.points.erase(vector<point_t>::Iterator(&pointset.points[rp]));
			#endif
		}
		else
		{
			count=0;
		}
	}
	#endif
}
void reset_out_of_bounds(pointset_t &pointset)
{
	int l=pointset.points.size();
	int i=0;
	while(l--)
	{
		double r=0;
		double x;
		x=pointset.points[i].x[0];
		x*=x;
		r+=x;
		x=pointset.points[i].x[1];
		x*=x;
		r+=x;
		x=pointset.points[i].x[2];
		x*=x;
		r+=x;
		r=sqrt(r);
		if(r>point_max_radius)
		{
			initialize_point(pointset.points[i]);
		}
		i++;
	}
}
void initialize_pointset(pointset_t &pointset)
{
	pointset.points=vector<point_t>(pointset_size);
	int num_points=pointset.points.size();
	while(num_points--)
	{
		initialize_point(pointset.points[num_points]);
	}
}

void define_materials()
{
	CLONE(&materials[0],m);
	m[0].shininess=32;
	/* white light component coefficients */
	SET4(m[0].specular,0.875,0.262,0.55,1);
	SET4(m[0].ambient,0.1,0.1,0.1,1);
	SET4(m[0].diffuse,1,.5,.5,1);
	SET4(m[0].emission,0,0,0,1);
	m[0].sides=1;
	m[0].orientation=1;
	m[0].rasterization_mode=material_t::LINE;
	m[0].interpolation_model=material_t::LINEAR;
}

void define_lights()
{
	/* ambient lighting for the entire scene */
	SET4(ambient_light,0,0,0,1);

	SET4(lights[0].x,10,10,10,0);
	SET4(lights[0].specular,1,1,1,1);
	SET4(lights[0].diffuse,1,1,1,1);
}

void reset_some(void)
{
	initialize_pointset(pointset);
}

void reset_all(void)
{
	draw_stereo=0;
	stereo_sep=0.1;
	alpha_magnitude=0;
	limd=1e-1;
	perturb=1e-4;
	do_perturb=1;
	enable_lights=0;
	define_lights();
	define_materials();
	draw_vector_field=0;
	draw_pointset_flag=1;
	animate_pointset_flag=0;
	vector_field_scale=0.1;
	define_vector_field();
	const_c=16;
	define_surface();
	draw_surface=0;
	draw_point=0;
	plane_scale=30;
	bool_draw_plane=0;
	plane_density=512;
	draw_param_curve=1;
	draw_tan_curve=0;
	draw_norm_curve=0;
	ortho=0;
	orscale=1;
	initialize_curves();
	calc_t();
	scale=1;
	scalex=1;
	scaley=1;
	scalez=1;
	mrotx=0;
	mroty=0;
	mrotz=0;
	mrd=3;
	eyex=5;
	eyey=5;
	eyez=5;
	cenx=0;
	ceny=0;
	cenz=0;
	upx=0;
	upy=0;
	upz=1;
	fov=30;
	aspect=1;
	znear=.01;
	zfar=1000;
	reset_some();
}

void cycle_rasterization_mode(material_t &mat)
{
	switch(mat.rasterization_mode)
	{
	    case material_t::POINT:
		mat.rasterization_mode=material_t::LINE;
		break;
	    case material_t::LINE:
		mat.rasterization_mode=material_t::FILL;
		break;
	    case material_t::FILL:
		mat.rasterization_mode=material_t::POINT;
		break;
	}
}

/* ARGSUSED1 */
void keyboard(unsigned char key, int x, int y)
{
	int i,n;
	switch(key)
	{
		case 'B':
			stereo_sep-=0.01;
			printf("%f\n",stereo_sep);
			break;
		case 'b':
			stereo_sep+=0.01;
			printf("%f\n",stereo_sep);
			break;
		case ' ':
			draw_stereo^=1;
			break;
		case 'w':
			cycle_rasterization_mode(materials[0]);
			break;
		case '1':
			/* view of xy-plane */
			eyex=0;
			eyey=0;
			eyez=5;
			upx=0;
			upy=1;
			upz=0;
			break;
		case '2':
			/* view of yz-plane */
			eyex=5;
			eyey=0;
			eyez=0;
			upx=0;
			upy=0;
			upz=1;
			break;
		case '3':
			/* view of xz-plane */
			eyex=0;
			eyey=5;
			eyez=0;
			upx=0;
			upy=0;
			upz=1;
			break;
		case '4':
			eyex=5;
			eyey=5;
			eyez=5;
			upx=0;
			upy=0;
			upz=1;
			break;
		case '5':
			vecfields[0].on^=1;
			break;
		case '6':
			vecfields[1].on^=1;
			break;
		case 'o':
			ortho^=1;
			break;
		case '':
			enable_lights^=1;
			/*
			DBUG(materials[0].diffuse[0]);
			DBUG(materials[0].diffuse[1]);
			DBUG(materials[0].diffuse[2]);
			DBUG(materials[0].diffuse[3]);
			*/
			printf("light:%f %f %f %f\n",
			       lights[0].diffuse[0],
			       lights[0].diffuse[1],
			       lights[0].diffuse[2],
			       lights[0].diffuse[3]);
			printf("diffuse:%f %f %f %f\n",
			       materials[0].diffuse[0],
			       materials[0].diffuse[1],
			       materials[0].diffuse[2],
			       materials[0].diffuse[3]);
			break;
		case 'D':
			universe_dampen-=0.001;
			break;
		case 'd':
			universe_dampen+=0.001;
			break;
		/*
		case 'd':
			scalez*=1.1;
			break;
		case 'D':
			scalez*=0.9;
			break;
		*/
		case 's':
			scale*=1.1;
			break;
		case 'S':
			scale*=0.9;
			break;
		case 'F':
			fov+=1;
			orscale*=1.1;
			break;
		case 'f':
			fov-=1;
			orscale*=0.9;
			break;
		case 27:
			exit(0);
			break;
		case 'x':
			mrotx+=mrd;
			break;
		case 'X':
			mrotx-=mrd;
			break;
		case '':
			mrotx=0;
			break;
		case 'y':
			mroty+=mrd;
			break;
		case 'Y':
			mroty-=mrd;
			break;
		case '':
			mroty=0;
			break;
		case 'z':
			mrotz+=mrd;
			break;
		case 'Z':
			mrotz-=mrd;
			break;
		case '':
			mrotz=0;
			break;
		case '':
			reset_some();
			break;
		case 'r':
			if(1)
			{
				int n=ELEMSOF(vecfields);
				while(n--)
				{
					int v=VFD;
					while(v--)
					{
						vecfields[n].n[v]++;
					}
				}
			}
			break;
		case 'R':
			if(1)
			{
				int n=ELEMSOF(vecfields);
				while(n--)
				{
					int v=VFD;
					while(v--)
					{
						vecfields[n].n[v]-=vecfields[n].n[v]>1;
					}
				}
			}
			break;
		/*
		case 'e': vecfields[0].n[0]++; break;
		case 'E': vecfields[0].n[0]-=vecfields[0].n[0]>1; break;
		case 'r': vecfields[0].n[1]++; break;
		case 'R': vecfields[0].n[1]-=vecfields[0].n[1]>1; break;
		case 't': vecfields[0].n[2]++; break;
		case 'T': vecfields[0].n[2]-=vecfields[0].n[2]>1; break;
	    */

		case 'g':
			if(0)
			{
				grav_threshold+=.01;
			}
			else
			{
				grav_mass*=2;
				fprintf(stderr,"grav mass %f\n",grav_mass);
			}
			break;
		case 'G':
			if(0)
			{
				grav_threshold-=.01;
				if(grav_threshold<0)
				{
					grav_threshold=0;
				}
			}
			else
			{
				grav_mass*=0.5;
				fprintf(stderr,"grav mass %f\n",grav_mass);
			}
			break;

		case 't':
		case 'T':
			if(1)
			{
				int f=1;
				if(key=='t')
				{
					f=-1;
				}
				int n=ELEMSOF(vecfields);
				while(n--)
				{
					int v=VFD;
					while(v--)
					{
						vecfields[n].lim[v][0]+=f*limd;
						vecfields[n].lim[v][1]-=f*limd;
					}
				}
			}
			break;
	/*
		case 'h': vecfields[0].lim[0][0]+=limd; break;
		case 'H': vecfields[0].lim[0][0]-=limd; break;
		case 'y': vecfields[0].lim[0][1]+=limd; break;
		case 'Y': vecfields[0].lim[0][1]-=limd; break;
		case 'j': vecfields[0].lim[1][0]+=limd; break;
		case 'J': vecfields[0].lim[1][0]-=limd; break;
		case 'u': vecfields[0].lim[1][1]+=limd; break;
		case 'U': vecfields[0].lim[1][1]-=limd; break;
		case 'k': vecfields[0].lim[2][0]+=limd; break;
		case 'K': vecfields[0].lim[2][0]-=limd; break;
		case 'i': vecfields[0].lim[2][1]+=limd; break;
		case 'I': vecfields[0].lim[2][1]-=limd; break;

		case 'n': vecfields[0].lim[0][0]+=limd; vecfields[0].lim[0][1]+=limd; break;
		case 'N': vecfields[0].lim[0][0]-=limd; vecfields[0].lim[0][1]-=limd; break;
		case 'm': vecfields[0].lim[1][0]+=limd; vecfields[0].lim[1][1]+=limd; break;
		case 'M': vecfields[0].lim[1][0]-=limd; vecfields[0].lim[1][1]-=limd; break;
		case ',': vecfields[0].lim[2][0]+=limd; vecfields[0].lim[2][1]+=limd; break;
		case '<': vecfields[0].lim[2][0]-=limd; vecfields[0].lim[2][1]-=limd; break;
	    */
		case ',':
			mag_resolution*=2;
			printf("mag_resolution=%d\n",mag_resolution);
			initialize_curves();
			break;
		case '<':
			if(mag_resolution>1)
			{
				mag_resolution/=2;
			}
			printf("mag_resolution=%d\n",mag_resolution);
			initialize_curves();
			break;
		case 'm':
			mag_current*=2;
			if(1)
			{
				fprintf(stderr,"mag current changed to %f\n",mag_current);
			}
			break;
		case 'M':
			mag_current*=0.5;
			if(1)
			{
				fprintf(stderr,"mag current changed to %f\n",mag_current);
			}
			break;
		case 'n':
			if(1)
			{
				add_points(pointset,256);
			}
			else
			{
				pointset_size<<=1;
			}
			break;
		case 'N':
			if(1)
			{
				remove_points(pointset,256);
			}
			else
			{
				pointset_size>>=1;
			}
			break;
		case '':
			vector_field_scale*=-1;
			break;
		case '':
			animate_curve^=1;
			break;

		case 'q': vecfields[0].normal_alpha_magnitude*=0.9; break;
		case 'Q': vecfields[0].normal_alpha_magnitude*=1.1; break;
		case 'a': vecfields[0].color[3]*=1.1; break;
		case 'A': vecfields[0].color[3]*=0.9; break;
		case '': alpha_magnitude^=1; break;


		case 0x0d:
			animate_pointset_flag=1;
			break;
		case '':
			draw_vector_field^=1;
			break;
		case '':
			draw_pointset_flag^=1;
			break;
		case '':
			do_perturb^=1;
			if(do_perturb)
			{
				printf("perturb on\n");
			}
			else
			{
				printf("perturb off\n");
			}
			break;
		case 'p':
			perturb*=1.1;
			break;
		case 'P':
			perturb*=0.9;
			break;
		case 'v':
			vector_field_scale*=1.2;
			break;
		case 'V':
			vector_field_scale*=0.8;
			break;
		case '':
			draw_param_curve^=1;
			break;
		case 0x13:
			/*
			draw_surface^=1;
			*/
			stop_pointset(pointset);
			break;
		default:
			printf("unbound key:0x%x\n",key);
			break;
		case '+':
			const_c+=1;
			printf("c=%f\n",const_c);
			break;
		case '-':
			const_c-=1;
			printf("c=%f\n",const_c);
			break;
	}

	glutPostRedisplay();
}


int main(int argc, char** argv)
{
	int window_width=512;
	int window_height=512;
	int fullscreen=1;
	int do_run_program=1;
	int mag_current_ma=(int)(mag_current*1000);

	struct parameter_set
	{
		char const *name;
		int *value;
		int *anti_value;
	};

	struct parameter_set parameter_set[]=
	{
		{
			"window_width",
			&window_width,
			&fullscreen,
		},
		{
			"mag_resolution",
			&mag_resolution,
			0,
		},
		{
			"mag_current",
			&mag_current_ma,
			0,
		},
		{
			"window_height",
			&window_height,
			&fullscreen,
		},
		{
			"fullscreen",
			&fullscreen,
			0,
		},
		{
			"pointset_size",
			&pointset_size,
			0,
		},
	};

	if(1)
	{
		if(argc>1)
		{
			int do_show_syntax=0;
			int i=1;
			int n=argc-i;
			while(n--)
			{
				int match=0;
				int k=sizeof parameter_set/sizeof parameter_set[0];
				while(k--)
				{
					if(strcmp(argv[i],parameter_set[k].name)==0)
					{
						if(n)
						{
							char *endptr=0;
							int value=strtoul(argv[i+1],&endptr,0);
							if(*endptr==0)
							{
								*parameter_set[k].value=value;
								if(1)
								{
									fprintf(stderr,"set value of %s to %d\n",parameter_set[k].name,*parameter_set[k].value);
								}
								if(parameter_set[k].anti_value)
								{
									*parameter_set[k].anti_value=0;
								}
								if(1)
								{
									match=1;
									k=0;
									n-=1;
									i+=1;
								}
							}
							else
							{
								fprintf(stderr,"invalid numerical value %s %s\n",argv[i],argv[i+1]);
							}
						}
						else
						{
							fprintf(stderr,"name argument requires value argument\n");
						}
					}
				}
				if(match)
				{
					if(0)
					{
						fprintf(stderr,"argument match: %s\n",argv[i]);
					}
				}
				else
				{
					fprintf(stderr,"argument mismatch: %s\n",argv[i]);
					do_show_syntax=1;
				}
				i+=1;
			}
			if(do_show_syntax)
			{
				int i=0;
				int n=ELEMSOF(parameter_set);
				fprintf(stdout,"%s program arguments:\n",argv[0]);
				while(n--)
				{
					fprintf(stderr,"\t%s <n>\n",parameter_set[i].name);
					i+=1;
				}
				do_run_program=0;
			}
		}
	}

	if(1)
	{
		mag_current=(double)mag_current_ma/1000;
	}

	if(do_run_program)
	{
		seed_random();
		reset_all();

		glutInit(&argc,argv);
		glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB);
		glutInitWindowSize(window_width,window_height); 
		glutInitWindowPosition (0, 0);
		glutCreateWindow(argv[0]);
		if(fullscreen)
		{
			glutFullScreen();
		}
		init();
		glutDisplayFunc(display); 
		glutReshapeFunc(reshape);
		glutKeyboardFunc(keyboard);
		display();
		glutMainLoop();
	}
	return 0;
}



